Diversity indicators using OBIS data (EEZ)
This notebook builds on the main notebook to generate indicator maps for EEZ only.
Create a discrete global grid
Create an ISEA discrete global grid of resolution 9 using the dggridR package and assign cell numbers to the occurrence data:
Calculate metrics
This calculates stations per 1000 square km. Make sure to adjust if the resolution is changed.
metrics <- occ %>%
group_by(cell, decimallongitude, decimallatitude) %>%
summarize(records = sum(records)) %>%
group_by(cell) %>%
summarize(records = sum(records), stations = n()) %>%
mutate(stations = stations / 2591.4 * 1000)
metrics## # A tibble: 122,014 x 3
## cell records stations
## * <dbl> <int> <dbl>
## 1 1 296127 510.
## 2 2 917 35.1
## 3 3 605 29.7
## 4 4 210 40.5
## 5 5 432 22.8
## 6 6 10 2.70
## 7 7 128 3.47
## 8 8 68 19.7
## 9 9 298 37.0
## 10 10 881 89.1
## # … with 122,004 more rows
Add cell geometries to the metrics table:
library(sf)
grid <- dgearthgrid(dggs, frame = FALSE, wrapcells = FALSE)
grid$cell <- names(grid)
grid_sf <- grid %>%
st_as_sf() %>%
st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=230"))
grid_cropped <- grid_sf %>% st_crop(c(xmin = -180, ymin = -85, xmax = 180, ymax = 85))
metrics <- merge(metrics, grid_sf, by.x = "cell", by.y = "cell", all.y = TRUE) %>%
#filter(st_intersects(geometry, st_as_sfc("SRID=4326;POLYGON((-180 85,180 85,180 -85,-180 -85,-180 85))"), sparse = FALSE)) %>%
st_sf()Load high seas shapefile
Use the high seas shapefile from https://www.marineregions.org/ to mask the high seas:
Create maps
library(rnaturalearth)
library(rnaturalearthdata)
library(viridis)
library(ggplot2)
world <- ne_countries(scale = "medium", returnclass = "sf")
ggplot() +
geom_sf(data = metrics, aes_string(fill = "stations", color = "stations"), lwd = 0.1) +
scale_fill_viridis(option = "inferno", na.value = "#000000", name = expression("OBIS: number of stations per 1000"~km^2), trans = "log10") +
scale_color_viridis(option = "inferno", na.value = "#000000", name = expression("OBIS: number of stations per 1000"~km^2), trans = "log10") +
geom_sf(data = world, fill = "#dddddd", color = NA) +
geom_sf(data = highseas, fill = "#ffffff", color = NA, alpha = 0.5) +
theme(
panel.grid.major.x = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.background = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position = "top"
) +
guides(fill = guide_colourbar(barwidth = 27, barheight = NULL, title.position = "bottom", title.hjust = 0.5)) +
xlab("") +
ylab("") +
coord_sf(crs = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")ggplot() +
geom_sf(data = metrics, aes_string(fill = "stations", color = "stations"), lwd = 0.1) +
scale_fill_distiller(palette = "Spectral", na.value = "#3288bd", name = expression("OBIS: number of stations per 1000"~km^2), trans = "log10") +
scale_color_distiller(palette = "Spectral", na.value = "#3288bd", name = expression("OBIS: number of stations per 1000"~km^2), trans = "log10") +
geom_sf(data = world, fill = "#dddddd", color = NA) +
geom_sf(data = highseas, fill = "#ffffff", color = NA, alpha = 0.5) +
theme(
panel.grid.major.x = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.background = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position = "top"
) +
guides(fill = guide_colourbar(barwidth = 27, barheight = NULL, title.position = "bottom", title.hjust = 0.5)) +
xlab("") +
ylab("") +
coord_sf(crs = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")ggplot() +
geom_sf(data = metrics, aes_string(fill = "stations", color = "stations"), lwd = 0.1) +
scale_fill_distiller(palette = "Spectral", na.value = "#3288bd", name = expression("OBIS: number of stations per 1000"~km^2), trans = "log10") +
scale_color_distiller(palette = "Spectral", na.value = "#3288bd", name = expression("OBIS: number of stations per 1000"~km^2), trans = "log10") +
geom_sf(data = world, fill = "#dddddd", color = "#000000", lwd = 0.2) +
geom_sf(data = highseas, fill = "#ffffff", color = "#000000", alpha = 0.5, lwd = 0.2) +
theme(
panel.grid.major.x = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.background = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position = "top"
) +
guides(fill = guide_colourbar(barwidth = 27, barheight = NULL, title.position = "bottom", title.hjust = 0.5)) +
xlab("") +
ylab("") +
coord_sf(crs = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")